Significant increase in graupel and lightning occurrence in a warmer climate simulated by prognostic graupel parameterization

There is little consensus among global climate models (CGMs) regarding the response of lightning flash rates to past and future climate change, largely due to graupel not being included in models. Here a two-moment prognostic graupel scheme was incorporated into the MIROC6 GCM and applied in three experiments involving pre-industrial aerosol, present-day, and future warming simulations. The new microphysics scheme performed well in reproducing global distributions of graupel, convective available potential energy, and lightning flash rate against satellite retrievals and reanalysis datasets. The global mean lightning rate increased by 7.1% from the pre-industrial period to the present day, which was attributed to increased graupel occurrence. The impact of future warming on lightning activity was more evident, with the rate increasing by 18.4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%\,\textrm{K}^{-1}$$\end{document}%K-1 through synergistic contributions of destabilization and increased graupel. In the Arctic, the lightning rate depends strongly on the seasonality of graupel, emphasizing the need to incorporate graupel into GCMs for more accurate climate prediction.

The observed and simulated global lightning distributions are also shown in Fig. 2, in which the widely reported land-ocean contrast (Fig. 2c,d) is evident, due to lightning being triggered not only by convection associated with atmospheric instability but also to frozen hydrometeors within clouds affecting charge separation 35 .Thus, the spatial distribution of lightning is not determined by CAPE alone (Fig. 2e,f).
The present study used the lightning scheme parameterized by both CAPE and precipitating ice hydrometeors coupled to prognostic graupel [Eq.( 3)], and the simulated spatial pattern of flash rate is similar to that of LIS/ OTD satellite observations (Fig. 2c,d).The model captured higher lightning flash rates well, particularly over Central Africa, the Amazon, the eastern coast of North America, India, and Southeast Asia.The flash rate is somewhat underestimated over the tropical convective zones, partly because of the coarse vertical resolution and because GCMs cannot explicitly resolve convective clouds.These issues may be addressed in a future study.
The global mean lightning rate increased by 7.1% from the PI to PD simulations (Table 1), due mainly to increased ice hydrometeors over the tropics (Supplementary Fig. 2), implying a signal of 'cloud invigoration hypothesis' 36 linking changes in lightning, graupel, and aerosols 14,[37][38][39] .The cloud liquid water path (CLWP), cloud ice water path (CIWP), and GWP increase in response to the change in aerosol from PI to PD, but the graupel is the most sensitive to this change (Table 1).This is likely because the increased cloud liquid and ice water contribute directly to the source of graupel.Given that CAPE is almost constant (Supplementary Table 1), the increased lightning rate is fundamentally attributed to the increased occurrence of graupel.The effective radiative forcing due to aerosol-cloud interactions ( ERF aci ; see "Methods" section) is − 0.62 W m −2 in the new MIROC6-CHIMERRA (Cloud-Hydrometeors Interactive Module with Explicit Rain and Radiation; see "Methods" section) model with prognostic graupel, ∼ 22% lower than that of the previous model without graupel (− 0.79 W m −2 ).This may be explained by the 'snow-induced ACI buffering hypothesis' 34 , which emphasizes the important role of the riming of underlying super-cooled liquid clouds by falling snowflakes and graupel at mid latitudes, effectively removing cloud water 40 even though anthropogenic aerosols are increased from PI levels.

Future changes in graupel and lightning
The response of lightning flash rate to a unit increase in global surface air temperature over higher northern latitudes (above 55 • N) is shown in Fig. 3a.An increased flash rate is predicted for most regions, with a mean increase over the Arctic (66 • N-90 • N) of 12.2% K −1 .The spatial pattern is inherently linked to that of graupel (Fig. 3b) and CAPE (Fig. 3c) responses [Eq.( 3)].More specifically, the distribution of the lightning change is generally similar to the change in CAPE, with the coincident increase in graupel canceling out the reduced lightning frequency due to reduced CAPE (i.e., stabilization) around the North Pole.
Although graupel and CAPE are the main controls on the lightning distribution in the model, the responses of meteorology and large-scale circulation are also important factors in lightning prediction because our simulations suggest that increased graupel is a result of incremental changes in CLWP and CIWP with increasing temperature (Table 1).Graupel size may be related to the cloud and environmental regimes (Fig. 1c), and future responses of graupel size to warming will also be independent of graupel occurrence and CAPE (Fig. 3d and Supplementary Fig. 3).The change in hydrometeor size is directly related to fall speed and thus the strength of the hydrological cycle, which may in turn be linked to charge separation within clouds 41 .Although beyond the scope of this study, an improved physical understanding of lightning occurrence and more elaborate parameterization may help to reduce model uncertainty and variability 12,16 .
The global mean sensitivities of graupel occurrence and CAPE to a unit temperature increase are comparable, + 7.9% K −1 and + 9.7% K −1 , respectively (Table 1).Over polar regions, however, the contribution of change in graupel occurrence is more evident (Table 1, and Supplementary Figs. 3 and 4).As a result of the combined effects of changes in graupel and CAPE, the global mean response of lightning flash rate is + 18.4% K −1 , higher than previous estimates of −5 to + 12% K −113, 42 .
This increase in lightning flash rate may be attributed mainly to different model performances in representing cloud vertical structures 17,43 , different treatments of convective clouds that are not resolved in GCMs 44 , or different microphysics process representations 19 .Although the model of the present study included the treatment of graupel in large-scale condensation, a fundamental physical consideration of future lightning changes with regard to convective clouds may have been overlooked 45 .The present results indicate the need for a prognostic treatment of precipitation, including high-density ice crystallites, for both cumulus and large-scale condensation; such a GCM is not currently available.
Table 1.Global and regional mean climatology of relative changes in lightning flash rate and associated parameters from pre-industrial aerosol levels (PI, year 1850) to present-day (PD, year 2000), and from PD through future warming (SST + 4K).A single asterisk (*) indicates statistical significance at the 95-99% confidence level; a double asterisk (**) indicates statistical significance at the > 99% confidence level.a CLWP, Cloud Liquid Water Path.b CIWP, Cloud Ice Water Path.c GWP, Graupel Water Path.d Q frz : Total precipitating ice hydrometer [see Eq. ( 4)]. e CAPE, Convective Available Potential Energy.Models and experimental configurations vary among studies in terms of, for example, future scenarios, CO 2 and aerosol emissions, horizontal/vertical resolution, and microphysical complexity.Although the responses of atmospheric circulation, surface forcing, and consequent feedback also differ between atmosphere-only and coupled-ocean configurations, there is no detailed documentation.Further study of model perturbation of statistics should be undertaken within a single-model framework.

Change in seasonal variation under the warming
There is evident seasonal variation in lightning activity 46,47 , and the possible effects of future warming on such variation are of interest but poorly documented.In this context, particular attention should be paid to polar regions, as the distinct occurrence of day and night cycles throughout the year in such regions provides an excellent opportunity for study of seasonal changes 48 .
Monthly time series of simulated graupel, CAPE, and lightning rate over the Arctic and Antarctic are shown in Fig. 4. In the Arctic region, graupel occurrence is higher in summer, particularly during June-September, and is lower in winter (Fig. 4a).This seasonality seems to be synchronized with that of CAPE (Fig. 4b) and therefore the resultant lightning flash rate (Fig. 4c), which is consistent with observational studies 47,49 .
Of note, although differences in CAPE between the PD and SST + 4K experiments were found only in their absolute values (Fig. 4b), differences in graupel occurrence between these experiments were also found in their seasonality (Fig. 4a) by shifting the maximal value to later in summer.Considering the relatively large difference in lightning flash rate between experiments for the period September-November (Fig. 4c), even though the CAPE did not change significantly, lightning activity must be influenced mainly by changes in graupel occurrence.
Such seasonal variation in graupel is also evident in the Antarctic (Fig. 4d) during February-April, but lightning rarely occurs poleward of 65 • S 5,50 due to the very stable environmental conditions in this region (Fig. 4e,f).Nevertheless, the Antarctic may be exposed to increased lightning under global warming (+ 26.4% K −1 ).
The results indicate that a reliable representation of seasonality in lightning activity under a global warming scenario is difficult to achieve unless the model considers the link between graupel and lightning, thereby emphasizing the need for the inclusion of prognostic graupel in GCMs.Furthermore, lightning occurrence may also be controlled by aerosol loading and type 38,49 .An understanding of interactions between microphysical aspects and lightning morphology at a fundamental process level 39,51 will be crucial in future studies.

Discussion
The response of lightning flash rate to aerosol perturbations and future global warming was considered using the MIROC6-CHIMERRA GCM.Most GCMs struggle to simulate lightning activity because graupel is ignored-a longstanding hindrance in lightning evaluation.Microphysical parameterization in CHIMERRA was therefore updated here by the inclusion of prognostic graupel.
The new microphysics scheme provided graupel distributions more consistent with those of cloud-resolving models and GPM/DPR satellite retrievals, and lightning flash rates consistent with LIS/OTD satellite observations (Figs. 1, 2).The global mean lightning flash rate is insensitive to PI through PD aerosol perturbations (+ 7.1%), but the impact of future warming on lightning activity is more evident.The global mean lightning rate is predicted to increase by 18.4% K −1 , more than in earlier studies 13,16,42 .
The increased lightning rate in the present model is attributed to destabilization and more frequent graupel, with the latter predominating at higher latitudes (Fig. 3 and Supplementary Fig. 3; Table 1).In particular, the lightning rate in the Arctic is strongly dependent on the seasonality of graupel (Fig. 4).The month of graupel maximum occurrence is delayed to late summer in response to future warming, resulting in a larger difference between PD and future climate during the autumn season (Fig. 4c).Changes in lightning rate in different climates associated with this seasonality would not be captured by a simplified lightning scheme depending only on cloud-top height, thereby emphasizing the need for improved ice microphysics in GCMs.Although some GCMs have included prognostic precipitation for large-scale condensation 21,22,52 , cumulus precipitation is still treated as diagnostically, and graupel and hail are generally not represented in GCMs.Future work on improving lightning projections will also lead to an improved understanding of lightning-produced NO x 53 and lightning- ignited wildfires 9 .
Variations among models are likely attributable to model performance in simulating cloud vertical structures, different treatments of convective clouds, and different microphysics process representations [16][17][18]54 . Theefore, model and scheme intercomparisons should be undertaken to elucidate the sources of uncertainty at the fundamental process level, with fixed model configurations.The present study may provide a benchmark for further experiments with different model configurations (e.g., involving CO 2 and aerosol emission scenarios, different horizontal/vertical resolutions, and atmosphere-only or ocean-coupled scenarios), microphysics, and lightning schemes.Overall, the degree of confidence in simulating future responses of lightning to global warming depends on prognostic graupel parameterization being included in microphysics-lightning interactions.

MIROC6 with CHIMERRA prognostic graupel parameterization
The latest version of the MIROC6 global aerosol-climate model was used.This differs from the same version in the Coupled Model Intercomparison Project Phase 6 (CMIP6) 55 in that it incorporates prognostic precipitation 22 .The model includes a two-moment large-scale microphysics scheme, the Cloud-Hydrometeors Interactive Module with Explicit Rain and Radiation (CHIMERRA), which considers both the mass and number mixing ratios of rain and snow.Precipitating hydrometeors also impact the radiation budget, so precipitation was radiatively active in the present model 34 .However, high-density ice hydrometeors (graupel and hail) are still not prognostic variables, which has been an obstacle to investigations of local orographic precipitation 29 and the link between graupel and lightning in GCMs.
Therefore, another version of the model was devised with prognostic graupel parameterization.The bulk graupel class was treated as a two-moment scheme, based on the following prognostic equations: where q g and N g are the mass and number mixing ratios of graupel, respectively; ρ a is air density ( kg m −3 ); u is the wind vector representing horizontal advection ( m s −1 ); v q g and v N are mass-and number-weighted fall velocities ( m s −1 ), respectively, representing vertical sedimentation; and S q g and S N g are respectively the source and sink tendencies ( s −1 ) of graupel with respect to deposition and sublimation, collection among other hydrometeors (i.e., accretion and riming), and melting, based on Reisner et al. 56 .Graupel was radiatively inactive.
The fall speed of larger graupel/hail may be > 10 m s −1 , which is computationally expensive during simula- tions of vertical sedimentation; the Eulerian solver requires short time-steps consistent with Courant-Friedrichs-Lewy (CFL) criteria.Therefore, the vertical sedimentation scheme was updated to a time-implicit precipitation sedimentation scheme 52 .The model was solved for microphysics with 60 s iterations 22 , allowing the sedimentation process to be solved at the same time-step, reducing computational costs.

Lightning parameterization and satellite data
The new parameterization with prognostic graupel enabled an evaluation of the frequency of lightning activity on a global scale.Among the various types of lightning scheme, it was convenient to apply one that explicitly considers graupel variation in considering the link between lightning rate and graupel morphology.
A scheme was chosen that combines frozen ice hydrometeors in the air column, Q frz ( kg m −2 ), and CAPE ( J kg −1 ) to parameterize the lightning flash rate, flight ( fl. m −2 s −1 ), following He et al. 23 : where α is a constant, and set to 2.67 × 10 −16 and 1.68 × 10 −17 for land and ocean, respectively, as applied by He et al. 23 for best performance among the various schemes; and Q frz is the column integrated precipitating ice hydrometeors (cloud ice, snow, and graupel) between the 0 • C ( z 0 ) and −25 • C ( z −25 ) isotherm levels: (1) where q i , q s , and q g are the mass mixing ratios ( kg kg −1 ) for cloud ice, snow, and graupel, respectively, which were prognosed in the model; and Q frz is a proxy for the charging rate due to collisions of graupel with other hydrometeors within clouds 57 .The parameterized lightning flash rate is highly variable, depending on CAPE and Q frz (Supplementary Fig. 5).Although Q frz includes cloud ice and snow in addition to graupel, the simulated lightning rate shows a good correlation with GWP (Supplementary Fig. 6).
For model validation, the Lightning Imaging Sensor and the Optical Transient Detector (LIS/OTD) are widely used 58,59 .This satellite dataset provides flash number density ( fl. km −2 ) with a 2.5 • × 2.5 • spatial resolution, and was used here in developing a climatology for the available period of July 1995-February 2014.

Experimental setup and evaluations
The following three sets of experiments were undertaken to investigate the sensitivity of lightning activity to aerosol perturbations and future warming: (1) a present-day (PD, year 2000) control experiment; (2) a pre-industrial (PI, year 1850) aerosol experiment; and (3) a future-warming experiment with a uniform 4 K increase in sea surface temperature (SST + 4K).Model resolution was 1.4 • × 1.4 • with 40 vertical levels, and an atmosphere- only configuration with year 2000 boundary conditions and prescribed climatological SST and ice.The model time-step was 12 min.A suite of simulations was undertaken for 11 years, and the subsequent 10 years used for analysis.
The impact of graupel treatment on the aerosol-cloud interaction (ACI) was assessed through ERF aci 60 quanti- fied in the PD and PI experiments as follows: where CRE is cloud radiative effect under clean-sky conditions 61 for the removal of contamination of aerosols (i.e., aerosol-radiation interactions or direct effects).
Future warming effects on graupel and lightning were evaluated by differences between experiments (1) and (3).Cloud, precipitation, and meteorological representations are highly model-dependent, so here the percentage change is reported in response to a unit surface air temperature perturbation ( % K −1 ), as reported in previous studies 13,42 .

Figure 3 .
Figure 3. Spatial distributions (above 55 • N) of future changes relative to PD ( % K −1 ) in (a) lightning flash rate, (b) graupel occurrence, (c) CAPE, and (d) graupel effective radius.Dotted regions indicate a statistically significant difference at the 95% confidence level.